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04 ' Practical methods to compute dipole strengths for a three-body system by using a dis- 

cretized continuum are analyzed. New techniques involving Green's function are developed, 
^ , either by correcting the tail of the approximate wave function in a direct calculation of 

. the strength function or by using a solution of a driven Schrodinger equation in a summed 

expression of the strength. They are compared with the complex scaling method and the 
Lorentz integral transform, also making use of a discretized continuum. Numerical tests are 
performed with a hyperscalar three-body potential in the hyperspherical-harmonics formal- 
ism. They show that the Lorentz integral transform method is less practical than the other 
methods because of a difficult inverse transform. These other methods provide in general 
comparable accuracies. 

1 Subject Index: 205 

, ^ , ' §1. Introduction 

\ Ah initio calculations of the strength function of a quantum system for a per- 

ly-^ ■ turbation play a very important role to understand its resonances and continuum 

Q^ , states. It is meant by ah initio that the calculation requires only a Hamiltonian and 

a physical condition of the system but is not based on specific model assumptions 
like a mean field model. This approach has attracted increasing interest with the 
advent of halo nuclei which are typical examples of a weakly bound system. These 
Q"^ I nuclei have only one or few bound states and most phenomena involving them are 

O ' properly understood by paying due attention to the role of the continuum and the 

resonances. 

Three-body continuum states are needed to calculate the energy spectra of the 
particles in the breakup of two-neutron halo nuclel^''^ as well as the three-body 
5^ \ partial decay widths of nuclei.'^' Though the construction of the three-body con- 

tinuum states is elegantly formulated in hyper spherical harmonics (HH) methods^ 
numerical computation generally involves heavy tasks and suffers from the problem 
of slow convergence .SI' ''3' 

Several methods have been proposed to resolve this problem within the HH 
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formalism P Another class of powerful approaches is to avoid a direct construction 
of the continuum but to utilize continuum-discretized states (CDS) as in the complex 
scaling method (CSM^ and the Lorentz integral transform (LIT^ method. Some 
of these methods are in fact applied to studying the breakup cross section and the 
electric dipole strength of the two-neutron halo nucleus ^He.'i^''^''^ 

The purpose of the present paper is twofold: The first is to propose a new 
method of generating three-body continuum in the HH formalism. The second is 
to assess the efficiency of this approach in a simple solvable problem by calculating 
the electric dipole strength function of a three-body system and comparing it to 
those of the CSM and LIT method. The basic vehicle of the present approach is 
a Green's function, the power of which is widely recognized in many problems and 
very recently proved in the phase-shift calculation of scattering between nuclei.^^ Ill 
tail behavior of the wave function obtained in a calculation using square-integrable 
basis functions is appropriately corrected with use of the Green's function. 

In Sect. [2] we briefly explain how the strength function is calculated not only 
in the CSM and LIT method but also in an approach based on the driven equation 
of motion (DEM). In Sect. [3] we present a Green's function approach to solve the 
three-body equation in the HH formalism. In solving both standard and driven 
equations of motion square-integrable bases are used in advantage owing to the 
Green's function. In Sect.[5]we treat a specific example of a three-body Hamiltonian 
which allows numerically exact calculations of the electric dipole strength function 
and thus enables us to compare to those obtained in various methods. Conclusion is 
drawn in Sect. O Appendix A presents a simple method of solving an inhomogeneous 
equation with an outgoing-wave boundary condition. The subtleties in the expansion 
of the complex scaled resolvent, which is at the heart of the CSM, are noted in 
Appendix B for the sake of completeness. The analytic form of the electric dipole 
strength for the three-particle continuum is discussed in Appendix C. 

§2. Methods of strength function calculation 

The physics we treat here involves the process that a system in its ground state is 
excited to three-body continuum states by a perturbation. Let H be the Hamiltonian 
of the system and W be the perturbation. We want to calculate the strength or 
response function given in lowest order perturbation theory as 

S{E) = ^ \{'F,\W\%)\^6{E, - E), (2-1) 

where is the ground state with energy Eq^ and is the excited state with 
energy E^: H]^q = Eq^q and H^^, = E^^y. The energy is measured from some 
standard value, e.g., a particle-decay threshold energy. When the excited state is in 
the continuum, the label u is continuous and the sum is replaced by an integration. 
In case where the energy Ey is degenerate, the sum also runs over the degenerate 
states. The continuum eigenstate ^y is normalized on an energy scale, i.e., {'I'y\^yi) = 

6{Ey -Ey,). 

It is convenient to derive an expression in which the summation in Eq. (|2-1|) 
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becomes implicit. This is possible if Eq. (|2-ip is rewritten as 

S{E) = {%\W^6{H - E)W\%), (2-2) 
which can be further rewritten as 

S{E) = --lm{%\W^g{E + ie)W\%). (2-3) 
vr 

Here 

GiE) = (2-4) 

is the resolvent, and use is made of the identity Q{E + ie) = FQ{E) — m5{E — H), 
where P stands for taking Cauchy's principal value of the integral. The outgoing- 
wave boundary condition is chosen by including ie in the resolvent, where e is a pos- 
itive infinitesimal value. It is noted that formulas for sum rule values, / E'^S{E)dE 
(n = 0, 1, ...), are easily derived with use of Eq. (|2-2p . 

The calculation of the strength function using Eq. ()2-3p was performed in Ref. ll4p 
to include the continuum effects of particle-hole excitations in the random-phase ap- 
proximation, and recently applied in the Hartree-Fock-Bogoliubov formalism.lI3''IS 
A new approach as well as some existing methods to calculate S{E) for more general 
cases are explained in the following subsections. 

2.1. Direct calculation with continuum discretized states 

The first and most obvious approach is to calculate expression ()2-ip directly. 
Because of the selection rules, only a limited number of wave functions must be 
computed at each energy. We shall use this technique to provide a reference cal- 
culation in Sect. 2] by computing the bound-state and scattering wave functions 
numerically with finite differences. 

However, our goal here is to show that the direct approach can be performed 
with square-integrable basis functions. Therefore, we can extend the Green's function 
technique developed in Ref. [TB|) to derive approximate scattering states directly us- 
able in Eq. ()2-ip . This direct method will be referred to as the continuum-discretized 
approximation (CDA). 

2.2. Driven equation of motion method 
Let ]F denote 

^ = g{E + ie)W'Fo. (2-5) 
The wave function is a solution of the driven Schrodinger equation 

{H - E)^ = -W% (2-6) 

with the outgoing-wave boundary condition in the asymptotic region. Thus we obtain 

S{E) = -lm.{^\W\^Q) = --lm.{^\H - E\^). (2-7) 
vr vr 

Solving Eq. (j2-6p with the outgoing-wave boundary condition is the main task in 
the DEM method. This approach is used to describe the double photoionization 
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of a two-electron atom!!^''!^ To cope with the boundary condition of the outgoing 
wave, the exterior complex scaling is often employed!^ Instead of using the exterior 
complex scaling, we will present in Sect. 13.21 and Appendix A a method of solving 
Eq. ()2-6p with use of the Green's function. 

2.3. Complex scaling method 

In the CSM the strength function is evaluated using the expression 

S{E) = --lm{%\W^U-\e)R{9)U{9)W\%), (2-8) 
vr 

which is readily obtained from Eq. (j2-3|) . where U{9) is an unbounded operator 
which transforms all the coordinates used to specify the system as x — > e^^x, and 
U-^{e) = {U{e))~^. Here R{9) is the complex scaled resolvent R{6) = U{9)g{E + 
ie)U-\9) = 1/{E - H{9) + ie) with H{9) = U{9)HU-^{9). 

Since any continuum eigenstate of H that has an outgoing wave in the asymptotic 
region is transformed into a function which damps at large distances within a suitable 
choice of 0, it is possible to expand the complex scaled eigenfunction over a set of 
linearly independent square-integrable basis functions <Pi{x)W^^^ 

H{9)^^{9) = E^{9)^^{9), (2-9) 



with 

^\9) = Y,CH0)Mx)- (2-10) 

i 

Here A is a label to characterize the eigenfunction. Substituting this expansion into 
Eq. (|2-9p enables one to obtain E^{9) and C^{9). The eigenfunction expansion of 
R{9) leads to the expression 

^ ' -K ^ E - E^{9) + ie' ^ ' 

A 

with 

V^{9) = {{'^^{9)y\W{9)\U{9)%), V^{9) = {{U{9)^Qy\W\9)\'^^{9)), (2-12) 

where W{9) = U{9)WU-^{9), W\9) = U{9)W^U-^{9), and U{9)^o is the solution 
of Eq. ()2-9p corresponding to the ground state, the eigenvalue E^{9) of which should 
be equal to Eq in principle. See Appendix B for the derivation of Eqs. ()2Tip and 
()2-12p . The accuracy of S{E) calculated with Eq. ()2-lip is tested by observing its 
stability against 9. See Ref. IHj) for the details, performances and references of the 
CSM. 

The contribution of a bound excited state to S{E) is not calculated from Eq. ()2Tip 
because its energy E^{9) becomes real in principle independently of 9, but is sepa- 
rately calculated from the original equation (|2-ip . 



Green's function method for strength function in three-body continuum 5 



2.4. Lorentz integral transform method 

Let us define the Lorentz transform of the strength function 

^(-) = r 7^^4^^dE = r --4^1— ,dE, (2.13) 



E^.. {E -z){E- z*) Je^^^^ {E - ErY + E] 

where z = Er + iEj is a complex energy and -Emin is a minimum energy from which 
S{E) begins to have strength. Substitution of Eq. ()2-ip leads to 

C{z) = l^'^-'^'fo^'' = {%\W^g{z*)g{z)W\%). (2-14) 

Thus £(z) reduces to the overlap, C,[z) = {^{z)\^{z)), of ^{z) = g{z)W%. The 
function ^{z) satisfies the following system of equations 

{H - z)^{z) = -W%. (2-15) 



This equation has the same structure as Eq. ()2-6p of the DEM method, but their 
solutions have quite different asymptotic behavior. Since C{z) takes a finite value in 
so far as z is complex, ^I'^z) has a finite norm, that is, it damps at large distances. 
Therefore ^{z) can be expanded in terms of square-integrable basis functions ^i{x) 
as 

^{Z) = Y,C^{Z)^^{X), (2-16) 

i 

and the coefficients Ci{z) are determined from the equation 

Y^mH - z\^j)C,{z) = -{<!>i\W\%). (2-17) 

An accurate calculation of C{z) is usually not a major problem, but determin- 
ing S{E) is the main task in the LIT method because it requires the inversion of 
Eq. (|2-13p . To perform the inversion, C{z) values are calculated for a number of Er 
values in a very wide interval for some chosen Ei value. This is necessary in order 
not to miss the sum rule of S{E), which is related to the integral of C{z) by 

/ S{E)dE = — C{z)dER. (2-18) 

One assumes a plausible form of S{E) which contains some parameters, and then 
determines those parameters so as to reproduce the C{z) data as accurately as possi- 
ble. One has to make sure that S{E) determined in this way is stable independently 
of the choice of Ej. See a review articl^ for the LIT method and the relevant 
literature. 

§3. Theory for three-body continuum 

3.1. Hyperspherical harmonics method 

We use a three-body system to compare the electric dipole strength functions 
calculated in various methods. To make this article self-contained, we give some 
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basic formulas which are needed to formulate the dynamics of three particles in the 
HH method. See Refs. [T]),[6]) for details. 

Let Aim, A2m, A^m and ^16,^26,^36 be the masses and charges of the three 
particles, where m is a nucleon mass and e is the electron charge magnitude. For 
the sake of simplicity we assume that two of the three particles are neutral, as in 
p + n + n and a + n + n, so that no Coulomb potential acts among the particles. 

Let H = T + V he the Hamiltonian of the three-body system with 

T = Ti+T2+n-T,^, 1^ = ^12 + ^23 + ^31+1423, (3-1) 

where Tj is the kinetic energy of particle i, and the kinetic energy of the center of mass 
motion, T^ui, is subtracted in the Hamiltonian. Here Vij is the nuclear interaction 
acting between ij pair and V123 is a nuclear three-body force. 
The Jacobi relative coordinates xi and X2 are defined as 

xi = ^/A^{rl-r2), X2 = 1/^12.3 ( "^^^^ ^ ^^^^ " , (3-2) 



where rj is the position vector of particle i, and Aij = AiAj/{Ai + Aj) and Aij^k = 
{Ai + Aj)Ak/{Ai + Aj + Afc) are the reduced mass factors. Let x^ denote the center 
of mass coordinate, x^ = (l/A) Yli=i ^ifi^ where ^ = ^1 + ^2 + ^3- 
The electric dipole operator of the system is given by 



= \j^'^Zie{ri - x-i)^ = ^J^idixi^ + ^22:2^) 

i=l 



(3-3) 



with 



dl = =(Zi^2 - Z2Ai), 

y^AiA2{Ai+ A2) 

^2 = " A {Zi + Z2)A, - ZMi + ^2)). (3-4) 

V V^l + ^2j^3^ 

Note that the quantity, + reduces to dl + dl = {Y,i{Zf/Ai) - Zif/A). 

The HH method is convenient to study the three-particle dynamics. The hyper- 
radius and hyperangle coordinates, p and q, are introduced as 



x\ + x\, a = tan "^(xi/x2). (3'5) 

Note that xi = psina, X2 = pcosa with < a < 7r/2. The five angle coordinates, a 
as well as 0i,(pi,92,4'2 of ;ri, X2, are denoted by il collectively. The volume element 
reads 

dridr2dr^ = j"^ p^dpdQdx^ (3-6) 

with 

/ A \ 
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where dxidx2 = p^dpdil, and d[2 = sin^ a cos"^ a da dxidx2- It is noted that the 
squared hyperradius is related to the mean squared radius, p'^ = Yli=i ~ 
x^)'^, or to the sum of the squared relative distances between the particles, p'^ = 

The kinetic energy operator T in Eq. (|3-ip reads 



2m \dx\ ~^ dx^) 2m ydp"^ ~^ pdp p^^ ^'^^) 



with the hypermomentum operator 

= -^-Acoi2a-^ + —^ll + —^Zl, (3-9) 
aa^ ua sm a cos^ a 

where ii,i2 are the angular momenta corresponding to the coordinates Xi,X2, re- 
spectively. 

The normalized eigenfunction of K? , called the HH, with eigenvalue K{K + 4) 
is given by 

^KLM.m = <t^%'\a)[Y,,{x^)Y,,{x2)]LM,, (3-10) 
where K is an integer called the hypermomentum, and 

4,^^^\a) = M^/^ sin^^ acos^^ aGn(li + ^2 + 2,£i + |; sin^ a), (3-11) 

where n is an integ er given by n — {^K — E-i — ^2)12 and the Jacobi polynomial G^ji is 
expressed in terms of the Gauss hyper geometric series as Gn{(i + ^2 + 2, ^1 + |; z"^) = 
F{-n, ei + £2 + n + 2,£i + ^;z^) and 



l 2{K + 2)r(£i + ^2 + n + 2)m + ^ + i) 

"V nim+n+Dim + Dr ■ ^'-'"^ 

The function -T^^lMl usually denoted y^j^^^ML literatures. Using the orthogo- 
nality of Gn leads to the following orthonormality relation 

/ i^K[\l^^)y^I^'hM'S^)d^ = Ke,Ki'M^K'5L,L'h4,Mi- (3-13) 

We also note the completeness relation 

6{Q -Q')= ^KLMSm^F'/^'uS^'W- (3-14) 

hl2KLML 
3.2. Solving with Green's function 

Ignoring the spin and isospin degrees of freedom of the particles, we focus on 
the spatial part of the wave function. The antisymmetry requirement on the wave 
function is also ignored. The inclusion of these causes no problem and will be detailed 
in a separate paper. 

We want to obtain the continuum state ^ for the equation of motion of type 

[H-E)^ = <P (3-15) 
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for a given Here E is a given energy from the three-particle threshold. There are 
two cases of our interest. In the first case, ^ = and \[^ is the continuum state to 
be used for 'F^, in Eq. (|2-ip . The second case is concerned with the DEM method, 
that is = — W^lZ'oj and we want to find such 'F that has the outgoing wave in the 
asymptotic region, as required in Eq. ()2-6p . In both cases we seek 'F with the angular 
momentum L and its projection Ml in an expansion in terms of the channel wave 
functions TcLMl = ^klmS^) 

^ = p-'^'Y.f-^P^^-LM,, (3-16) 

c 

where c = {£ii2K) stands for a set of channel labels. The aim is now to obtain the 
hyperradial function fdp)- 

Owing to the orthonormality (|3-13|) . fdp) is related to \F as 

fdp)=p'/'{J'cLM,m. (3-17) 

Substituting Eq (j3T6p to Eq. (|3-15p and projecting it to the channel c, we obtain 
^2 + ^ )fc{p) = j^{zc{p) + Cc{p)}, (3-18) 

with 

Zdp) = p'/\T,LMjVm, Up) = -p'/\TcLM,m, (3-19) 

where A;^ = 2mE/h'^. Note that the function z^p) is unknown. By using the 
expansion 

Up) = Y.(^-LM,\V\Tc'LM,)fAp) = Y.^cAp)fc'{p) (3-20) 

c' c' 

as usually done, Eq. ()3-18p becomes a set of coupled equations for fdp), but we keep 
the form (|3-18p in order to develop our approach. The function Qdp) known: It 
is either zero for ^ = or p^^"^ {J^cLMl\Wp\^q) for = —W^^q. It is important to 
realize that both of zdp) and Cc(p) are finite-ranged, that is, they vanish for large p. 
This is because V is finite-ranged and the ground state wave function Fq is spatially 
confined. In actual cases zdp) is known to decrease in power law (n > 3) even 
for a pairwise short-ranged nuclear interaction.'^ Moreover, when the Coulomb 
potential is present, the coupling between the different channels persists.^ These 
problems will make zdp) decrease very slowly. 

Let vdp) and hdp) denote, respectively, the regular and irregular solutions of 
the homogeneous equation with the right-hand side of Eq. (j3-18p set to zero. They 
are given in terms of the Bessel functions of the first and second kinds 

vdp) = {kpf'^JK+2{kp), hdp) = {kpf'^YK+2{kp). (3-21) 

They satisfy the Wronskian relation, W{vc-,hc) = vdp)h'dp) — Up)^c{p)= 2k/iT, 
where h'dp) = {d/ dp)hdp) etc. Now we discuss how to solve Eq. ()3T8p for the two 
cases. 
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(1) CD A case 

The function fc{p) must be regular at p = 0. By noting Cc{p) = 0, the formal 
solution of Eq. (|3-18|) can be written as 



2m 1'°° 

Act^c(p) + G,{p,p')z,{p')dp', (3-22) 



where Ac is a constant yet to be determined. The Green's function Gc is a solution 
of the equation 



dp'^ p^ 



2> +k2\G,{p,p') = 5{p-p\ (3-23) 

and it is given by 



Gc{p,p') = ^vc{p<)h,{p^), (3-24) 
where /0< (/0>) is the lesser (greater) of p and p'. The asymptotic behavior 

f!^^ip) AeK(/o) - tan 6Mp)} (3-25) 

p— ^-oo 

provides the phase shift 5c by 

vrm f°° 

tan 5c = - f^2j^x J '"c{p)zc{p)dp. (3-26) 

The normalization of f2^{p) will be discussed in Sect. 13.31 

We have to know Zc{p) and determine Ac. A basic idea for resolving this prob- 
lem is to make use of the CDS following Ref. The CDS are obtained by 
diagonalizing the Hamitonian in a certain square-integrable basis set which can ac- 
curately describe the wave function in the interaction region. As mentioned above, 
Zc{p) is finite-ranged, so that the ill tail behavior of causes no problem for 
evaluating Zc{p) reliably. We only need to require that is accurate in the in- 
ternal region where the potential V is effective. Thus we replace Zc{p) in Eq. (|3-22p 
by z^^{p) = p^/'^ {FcLML\yW^^) ■ Moreover fc{p) calculated from Eq. (|3-17p using 
denoted fc^{p) = p^^"^ {FcLMi^W'^^) ■, is expected to be already accurate in 
the internal region. Comparing it with f^^{p) in the internal region enables us to 
determine Ac- The procedure to determine Ac is as follows. Let [0, ad be the internal 
region. Choosing M sampling points {pi,p2, ■ ■ ■ , Pm) {Pi £ [Oj^c])) we determine Ac 
by the least squares fitting 

minimize over Ac : J][/c^^(pi) - f?^{pi)?. (3-27) 

i 

In this way we have the continuum state at hand and hence can calculate the strength 
function directly. 

As seen above, the energy E of Eq. ()3-15p cannot be chosen arbitrarily in the 
CDA, but has to be set to the discretized energies determined by the diagonalization. 
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(2) DEM case 

In the DEM method we have to solve Eq. (j3-18p with Cc{p) 7^ 0. A general 
solution that is regular at p = and has an outgoing wave in the asymptotic region 
reads 

fciP) = J^l Gt{p,p'){zc{p') + Up')}dp'. (3-28) 
The Green's function satisfying the outgoing-wave boundary condition reads 

Gt{p,p') = ^.vc{p<)ht{p>), (3-29) 

with 

ht{p) = vc{p) + ihc{p). (3-30) 

The function Zc{p) is again unknown. We note that fdp) and Zc{p) are mutually 
linked by Eqs. ()3-20p and ()3-28p . A method to determine fdp) is as follows. Because 
both Zcip) and Cc{p) are finite-ranged, that is, negligible for, say p > ac, it is found 
that fc{p) takes the form for p > Uc 

Up) = A,htip), (3-31) 

with a constant 

^c = ^J^ Mp){zc{p) + Up)] dp. (3-32) 

Therefore we need to determine fdp) in the region [0, ad in such a way that it joins 
Eq. ()3-3ip smoothly at p = ac- This can be performed in a basis expansion method 
as explained in Appendix A. In order for this method to work, we have to make sure 
that the strength function is stable with respect to the change of ac- 

3.3. Normalization of continuum states 

To discuss the normalization of the continuum state in CDA case, we first con- 
sider the plane wave (|27r)~^e**'i'^i^**'2'^^, which is expanded as follows^ 

(3-33) 

where k'^ = kf + /^l and fi^ denotes the five angles constructed from ki and ^2 hi 
exactly the same manner as f2. Using this expansion and the completeness rela- 
tion ()3T4p . we have 

pJK+2{kp)JK+2{k'p)dp = -5{k-k'). (3-34) 

The general form of the properly normalized free-wave that has energy E and the 
angular momentum L and its projection Ml is 

^ITm, =CkYl C'/-{k){kp)-'JK+2{kp)T'j^t'M,i^)^ (3-35) 
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where the amphtudes C^j^^^{k) satisfy the condition Yleie2K = 1 and Ck 

is a normahzation constant on the energy scale, {^k'^Mj^kLAiL^ ~ ^{E — E'). Using 
Eq. (fTOl) and 5{E - E') = 6{{h^/2m){e - k'^)) = {m/h'^k)6{k - k'), Ck becomes 



Ck = J''^k\ (3-36) 
n 

With reference to the above result, the normalization of f^^ may be chosen as 

cos 5r 



CD I 



A. 



/c^^(p), (3-37) 



where ||/p°|| = [Ji^ifc^{p))^dpY/^ is calculated from Eq. ([3071) using 0/^° and has 
the property X]cl|/^^|P=1 provided that tp'CD jg normalized. A continuum state 
^kLML normalized on the energy scale is given by 

^^Fm, = P-'^' E AT(p)-^cLM.. (3-38) 

c 

It should be stressed that the CDA method of constructing the continuum state 
does not require solving the coupled equations for fdp) but only needs the CDS. 
However, it does not allow chosing an arbitrary energy. Then the discretized state 
is expanded into the channel components and the tail behavior of each hyperradial 
part is readily corrected with the Green's function. 

§4. Specific examples 

The masses of three particles are set equal, j4i = = ^3 = 1, in unit of 
h? /m = 41.47106 Me Vfm-^. One of the particles has charge e and others are neutral. 
The three particles are assumed to interact via a hyperscalar potential which depends 
on the hyperradius only 

V = V{p). (4-1) 

As commented below Eq. (|3-7|) . p scales the size of the system, and the potential 
depending on p is considered a special three-body force. No channel coupling oc- 
curs for this potential in the HH formalism, and the three-body problem actually 
reduces to an easily solvable potential problem. The ground state consists of a sin- 
gle channel cq = (0,0,0), and its wave function takes the form /9~^/^/co(/5)-Fooo(^)- 
The excited states with L = 1 which are excited from the ground state by the elec- 
tric dipole operator have two channels, ci = (0,1,1) and C2 = (1,0,1), and they 
are degenerate in energy. Their wave functions are given by /9~^/^/ci (/9)J^|^^ (I?) 

and p~^/'^ fc2{p)J^liM^{^)y respectively. Note that the electric dipole operator p-3p 
acting on the ground state leads to 

UcAp) = UcM = -^Pfcoip). (4-2) 
"1 «2 vovr 

Thus the electric dipole strength function is easily evaluated if the hyperradial func- 
tions fc{p) are obtained. 
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Table I. The potential parameters and the ground state properties. Set 4 potential is three-ranged. 
Energy and length are given in units of MeV and fm. 



Set 


Vo 




Eo 




1 


-110 


0.16 


-17.6 


1.39 


2 


-90 


0.16 


-8.95 


1.57 


3 


-75 


0.16 


-3.49 


1.84 


4 


-610 


0.25 


-38.4 


0.938 




570 


0.16 








-200 


0.10 







The function fdp) is determined from Eq. ()3-18p . Since Zc{p) of Eq. ()3-19p 
reduces to V{p)fc{p) for V = V{p), Eq. ()3-18p simplifies to 

^2 W^^f"^ + ^ I ^'^P^ = 'W^^^P^- ^ 

In the CDA case of Cc{p) = we can solve this equation for both bound and contin- 
uum states with high precision using, e.g., the Numerov method. With this solutions 
we can numerically obtain the exact electric dipole strength function, and assess the 
various methods described in the previous section by comparing their strength func- 
tions with the exact one. In the DEM case where Cc(p) does not vanish, the above 
equation can easily be solved using the Bloch operator formalism^ 'i^'i^ or the 
method of Appendix A. 

For the CSM to be applicable, the potential V{p) should satisfy analyticity, 
and a square well potential or a Woods-Saxon potential with a small diffuseness 
parameter has to be avoided. The form factor of V{p) is assumed to be Gaussian, 
Vbexp(— Kp^). To generate the electric dipole strength functions of different shape, 
we chose four sets: Three of them are one-ranged Gaussian and the other is three- 
ranged. The strength and range parameters are listed in Table HI The ground state 
energies Eq and the root mean square value, ^^(/O^, of the ground state are also 
listed in the table. The Set 4 potential combined with the effective centrifugal barrier 
35fi^/{8mp^) {K = 1) has double minima at about p = 1.4 and 5.0 fm, producing a 
strength function of complex shape. 

We used the functions 

= /9^+^/'exp(-ia,/,2) (4-4) 

as square-integrable bases for an expansion of fdp)- The parameters were chosen 
in a geometric progression, Oj = l/(6oP*~^)^ with 6o = 0.1 fm and p = 1.3, to cover a 
wide p-space. The maximum number of terms was increased up to i = 31 to assure 
the convergence of the calculation, particularly in the LIT case. The electric dipole 
strength function shown below is not S{E) but s{E) = S{E)/{df + d^) nnit of 
fm^ MeV~^ 

4.1. Result with CDA 

Equation (j4-3p with Ccip) = is solved in a combination of basis functions (f){ai) 
in order to get the CDS, Figure [T] compares the phase shifts for K = 1{L = 1) 
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Fig. 1. (Color online) Comparison of the P-wave {K — 1) phase shifts calculated with Numerov 
and CDA methods for the potential sets of Table U] 



> 

CD 



E 



1 1 1 

Numerov 




CDA 




CDS . 

1 1 1 «_ 



> 



E 

iir 





15 



20 



15 



20 



Fig. 2. (Color online) Comparison of the electric dipole strength function between Numerov and 
CDA methods. Dotted vertical lines stand for the strength calculated with CDS. Potentials: 
Set 1 (Left upper), Set 2 (Left lower), Set 3 (Right upper). Set 4 (Right lower). 



14 



Y. Suzuki, W. Horiuchi, and D. Baye 



between the Numerov and CDA methods. The CDA phase shift shown in symbols 
is calculated from Eq. (j3-26p . None of the sets supports a bound state with L = 1, 
but Set 1 potential produces a sharp resonance around E = 2.5 MeV. The resonance 
is shifted to higher energy and becomes broader for the potentials of Sets 2 and 
3. The CDA calculation reproduces the exact phase shifts fairly well for all four 
potentials. This confirms that the wave function constructed from the discretized 
state with the help of the Green's function can describe the continuum accurately. 
The CDA can thus reproduce the strength functions without question, as shown in 
Fig. [21 Correcting the tail behavior of the discretized state is therefore very useful to 
predict the exact strength. The dotted vertical line in Fig. [2] stands for the strength 
calculated from the CDS, which ignores the continuum effect. The CDS strength 
is in reasonable correspondence to that of CDA near the sharp resonance region at 
about E = 2-3 MeV for Set 1 potential, but it generally exhibits a strong deviation 
at higher energies. 

4.2. Result with DEM 



To solve Eq. ()4-3p in the DEM method following Appendix A, we again use 
the same basis functions 4>{ai). Figure [3] displays the s{E) values for Set 4 and 
Set 2 potentials obtained in the DEM calculation. The strength for Set 4 becomes 
stable for Qc > 10 fm, while a larger ac is required for the Set 2 case to obtain the 
stability. The performance in DEM is very satisfactory for all the potential sets, that 
is, accurate strength functions are obtained independently of the potentials once Oc 
is chosen to be large enough. 



> 



E 



0.07 
0.06 
0.05 
0.04 
0.03 
0.02 
0.01 




T r 

Numerov 




10 
E [MeV] 



Fig. 3. (Color online) Comparison of the electric dipole strength function for Set 4 (Left) and Set 
2 (Right) potentials between Numerov and DEM calculations. 



4.3. Result with CSM 

In the CSM calculation p is transformed to e*^/9. The diagonalization of the 
rotated Hamiltonian is performed using the basis functions </)(aj). Figure S] displays 
the complex eigenvalues E^{6) for Set 4 potential for some 9 values. We observe 
two eigenvalues which are rather stable with respect to the change of 6. One of 
them corresponds to a sharp resonance at about E = 3 MeV in Fig. [H and another 
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to a very broad peak around 7 MeV. The latter shows up for 9 > 15°. Figure [5] 
displays how the strength function chang function of 9. The scaling angle 9 

is changed up to 25° to see the dependence of s{E) on 9. The CSM reproduces the 
exact strength very well when 9 is taken in the range 15° ^ 9 < 20° which covers the 
two stable eigenvalues noted above. Increasing 9 beyond 25° begins to deteriorate 
the agreement with the exact strength function. 



> 

0) 



" -3 



LU 



3 1 1 1 1 r 

+ + + 



+ 



o 



□ 



9=5° 
10° 
15° 
20° 
25° 



J I I L. 



1 2 3 4 5 6 
Re(E) [MeV] 



7 8 



Fig. 4. (Color online) Complex eigenvalues with K = 1, L = 1 for Set 4 potential as a function of 
the scaling angle 9. 
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Figure [6] plots the relative error, {sexa,ct{E) — s{E))/ Sexa.ctiE), of the strength 
functions calculated with the CSM and DEM models, where Sexact(-E') stands for 
the strength function calculated by the Numerov method. Both CSM and DEM 
reproduce the exact strength very well (within 1 %) over the wide energy range 
except that the CSM tends to give large errors near the threshold energy. The reason 
for this is understood from the fact that a correct energy dependence of s{E) for E 
close to zero is not manifestly guaranteed in the CSM. The DEM can incorporate 
such a behavior because it solves the driven equation of motion directly. 



0.02 



0.01 



.> 
tr 



-0.01 



-0.02 



10 
E [MeV] 



15 



20 




0.02 



0.01 



0.01 



0.02 



Fig. 6. (Color online) Relative error of the electric dipole strength function for Set 2 (Left) and Set 
4 (Right) potentials. 



4.4. Result with LIT 

The same type of bases as in CDA is used to calculate the Lorentz transforms 
C{z). Since we need C{z) in a wide Er region, the falloff parameters at in Eq. ()4-4p 
are taken in a region wide enough to describe various shapes of the function 1^(2;). 
To make the inversion from C{z) to S{E), it is convenient to express S{E) in terms 
of some plausible functions which the dipole strength of three particles is expected 
to take. As discussed in Appendix C, S{E) should show E"^ dependence near E = 0, 
so we assume the following form 

N 

s{E) = E-'Y,Cnexj>(^-^E), (4-5) 

n=l 

which is often employed in the literature.'^ The C{z) values calculated as the norm 
of 'I^{z) are fitted from those calculated from the above s{E) using Eq. ()2-13p . The 
coefficients Cn are determined by the least squares fitting. The parameter a and the 
maximum number of terms are varied to obtain "converging" s{E) as much as 
possible. Figure [3 compares the two C{z) for some Ej values. 

Figure [5] displays the corresponding s{E) functions. Though the fitting of C{z) 
appears to be very satisfactory except for the case of small Ej, the resulting s{E) 
functions differ from each other, exhibiting some deviations from the exact strength 
function in the energy region below the peak of the strength. Particularly the 
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Fig. 7. (Color online) Comparison of the Lorentz transform as a function of Er for Ej = 
5, 10, 15 MeV. Solid lines are obtained as the norm of ^'(z), while dotted lines are the inte- 
gral transform of Eq. H4-5P with the Lorentz kernel. Set 3 potential is employed. 



0.18 




5 10 15 



E [MeV] 

Fig. 8. (Color online) Comparison of the electric dipole strength function between the Numerov 
method and the LIT methods calculated for Ej = 5, 10, 15 MeV. Set 3 potential is employed. 



strength near the threshold shows an oscillatory behavior and becomes even neg- 
ative. The latter is unavoidable in general because the expression (|4-5p does not 
guarantee the positive definiteness. We tested another form 

/ N ' \^ 

s{E) = E'iS2C'^eKp[-^E)\ , (4-6) 

which is always positive. In this case determining the coefficients is not so easy. 
We used the amoeba routine—' to search for a minimum of a function. Since it is 
hard to increase N in this routine, no better result was attainable. 
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In such a case that the strength function has a very narrow peak or complex 
shape as a function of energy, the assumption of S{E) becomes in general much 
less trivial. Though there is some effort jS^'ESJiEZli no convincing form of the strength 
function or inverting method is available. 

§5. Conclusion 

The electric dipole strength functions provide interesting possibilities of compar- 
isons between theory and experiment. They require accurate theoretical calculations 
which become very difficult for a three-body system. In recent years, the three-body 
continuum has started to play a vital role in the study of the excitation mechanism of 
two-neutron halo nuclei. In the present paper we have developed and compared prac- 
tical methods to compute dipole strengths for a three-body system with a discretized 
continuum. Using discrete square-integrable states has the important advantage that 
powerful techniques developed for bound-state variational studies can be employed 
such as the versatile expansions in Gaussian states. 

The new techniques to compute the strength function with a discrete basis are a 
direct approach with discrete states (CDA method) and the calculation of a summed 
expression of the strength function involving a solution of an inhomogeneous driven 
equation of motion (DEM method). They both make use of Green's functions. 
We have compared them with the complex scaling method (GSM) and the Lorentz 
integral transform (LIT), also making use of a discretized continuum. An apparent 
advantage of the GSM and LIT is that strength functions can be obtained from 
states with unphysical asymptotic behaviors, i.e. asymptotically vanishing functions. 
However, the Green's function technique leads to physical asymptotics. In the CDA, 
the correct asymptotic behavior of scattering states is first obtained, allowing a direct 
calculation of the dipole strength. In the DEM, the solution of the driven equation 
of motion has a three-body outgoing-wave behavior. 

Numerical tests have been performed with an hyperscalar three-body potential 
in the hyper spherical-harmonics formalism. They have shown that the LIT method 
is not as mature as the other ones. The problems encountered with the LIT are 
not due to the discretized continuum but to the difficulty of accurately inverting the 
obtained transform. Any significant progress in this inversion would improve the 
present results. All other methods have provided comparable accurate results. The 
GSM presents however some accuracy problems at very low energies. 

The present approaches are promising tools for the study of realistic three-body 
strength functions with full account of channel couplings. In the three-body contin- 
uum, an infinity of open channels occur at each total energy. They correspond to the 
various ways of how this energy can be shared between the particles. The necessary 
extensions will also have to deal with Goulomb potentials for which couplings occur 
even at large distances. 
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Appendix A 

Solution of an inhomogeneous equation with a boundary condition 



The aim of this appendix is to solve Eq. (|3-28|) together with (jc{-2Up . where 
Zc{p) and Cc{p) are assumed to vanish for p > ac- From Eq. (|3-3ip . the logarithmic 
derivative of fdp) at p = Oc is a known constant given by 

„ acfcjac) achfjac) ^^ 
tc{ac) hj[ac) 

We thus need to determine fc{p) in the interval [0,ac], which we call /c"*(p), with 
the constraint that its logarithmic derivative at /? = Oc is Be- 

We try to obtain f^^^{p) in an expansion with some basis sets {(I)ci{p),4'c2ip), ■ ■ ■ , 

4>cn,{p)) 

ric 

fT\p)=Y.^-^c^{p) (A-2) 
1=1 

with 

l^ciip) = 4>ci{p) - lci4>co{p)- (A-3) 

Here (pcoip) is an auxiliary function introduced to satisfy the condition for the log- 
arithmic derivative at p = Oc- In fact the logarithmic derivative of fc^^{p) at p = ac 
becomes Be provided that jci is chosen as 

Bc<t)ci{ac) - ac(l)ci{ac) 
Bc(Pco{ac) - Oc0co(«cj 

The equation to determine the coefficients Xd reads 

J2 J2{(^cMc'j)^c,c' - &JMcc'\^c'j)}Xc'j = &JSc), (A-5) 

c' j=l 

where the round brackets indicate that the integration is to be done in [0, ac]: 



) = / {'Pciip)y'Pc'iip)dp, 
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_ _ 2m f"''^ l'°'<^' — — 

(</.,,|A4,,H'/'c',) = ^ j J [<l^,M*Gt{p,p')V,Ap')<l^c'j{p')dpdp', 

{^M = -^J^J^ {cP,,ip)rGtip,p')Up')dpdp'. (A-6) 

The amplitude Ac is determined from {ac) / {ac) ■ 

Appendix B 

Resolvent in complex scaling method 

The complex rotation U{9) in the CSM is defined by 

U{e)^{x) = exp Qi/e) ^xe''), (B-1) 

where / is the degree of freedom. (/ = 2 for a three-body system). Using U {9)U~^{9) 
= 1, it is easy to show that U~^{9) = U{—9). A key point is that this complex 
transformation makes the outgoing wave damp asymptotically. For a real potential 
V{x), the Hermitian conjugation of H{9) is H{9)^ = H{-9) = H{9)*, where * 
denotes the complex conjugation. Note that the eigenfunctions of H{9), ^^{9) and 

(9) of Eq. ()2-9p . with different labels A and A' are not orthogonal in general. 

Together with Eq. (j2-9p . we consider an accompanying eigenvalue problem 

H{9)${9) = E{9)^{9), (B-2) 

with H{9) = H{9)* = H{9)^ . The solution of this equation is labeled by the same A 
as that of Eq. (j2-9p . and we may choose the solution as follows: 

^^{9) = {^^{9))*, E^{9) = {E^{9)y. (B-3) 

We can show that both sets of {^'^{9)^ and {^'^{9)^ are biorthogonal,^'^' that is, 

{^^\9)\^\9)) = 6x^x., (B-4) 

if the normalization of ^■^{9) is chosen to satisfy 

j {W\9)fdx = 1. (B-5) 

To prove Eq. ()B-4p . we start from 

{§^'{9)\H{9)\^^{9)) =E^{9){^^'{9)\^^{9)). (B-6) 

The left-hand side is reduced to 

$r^'{9)\H{9)\^^{9)) = {H{9)^^'{9)\^^{9)) = E^' {9){^^' {9)\^^{9)). (B-7) 

Thus for E^{9) / E^' {9), {^^' {9)\^^{9)) = 0, which, together with Eq. (ITO . leads 
to the biorthogonality relation ()B-4p . 



Green's function method for strength function in three-body continuum 21 



It follows from the biorthogonality that the resolvent can be expanded as 

Substitution of this R{0) into Eq. (j2-8p leads to the strength function ()2-lip . but 
V^{e) and V^{e) are defined by 

v^{e) = {^^{e)\u{9)w^o), v^{e) = {^Q\w^u~^{e)^^{e)). (b-9) 

These are shown to be identical to those given in Eq. (j2T2p . The equality for 
T)^{6) is trivial from the definition of W{9). For the case of T>^{9), we only need 
to show that (^o\U~^{9) and (([/(0)!P'o)*| are identical within a ^-independent phase 
which may be chosen to be unity. This is justified by comparing the normalization 
condition, ((f7(0)!P'o)*|C/(0)lZ^o) = 1, with the trivial normalization condition for iFq, 

{^o\u~\e)u{e)\^^) = 1. 

Appendix C 

Analytic form of electric dipole strength function in three-body continuum 



The form of S{E) as a function of E is vital in inverting Eq. (j2-13p which 
relates S{E) to C{z). Its form may be discussed by examining the matrix element 
{\I^u\Wfi\'l^o) . To this end we assume that the continuum state is approximated 
by the free wave (|3-35p . First we consider a case where E is close to zero. Since 
K = 1 in the present case, the free wave has /c^-dependence for small k, that is, the 
matrix element behaves like oc E^^'^. Hence S{E) has an E^ dependence near the 
threshold. Note that this dependence is different from that of a two-body continuum, 
where the radial part of the free wave normalized on the energy scale is given for the 
partial wave i by 

^Mtr). (CI) 

Here fi is the reduced mass of the two particles. Thus the electric dipole matrix 
element for small k scales as /c^/^ oc E^^^ for the P-wave. 

To know a more general form of S{E), let us assume that the hyperradial part 
fcoip) of the ground state iP'o is approximated in terms of a combination of Gaussians, 
p^/^ exp(— ^ap^), or Exponentials, p^/^ exp(— 6p). The ^-dependence of the matrix 
element becomes /c^ exp(— A;^/2a) for the Gaussians, while it is given by k^{k'^/b'^ -\- 
2^-9/2 £qj, ^j^g Exponentials. The strength close to = is included in these 
functional forms. The analytic form of S{E) we are looking for is therefore 

5(i?) = ii;3('^C„e-"^) , (C-2) 
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with Q„ = m/{h^an) for the Gaussian case, or 

with = 2m/(/i^6^) for the Exponential case. 

Because the actual continuum deviates from the free wave through the inter- 
action among the particles, neither Eq. (|C-2p nor Eq. (|C-3p is exact. However, a 
suitable choice of the coefficients C„ may simulate S{E) accurately. The Lorentz 
transform is calculated from Eq. (jC-2p or Eq. (|C-3p according to Eq. (|2-13p . The 
coefficients C„ are then determined by least squares fitting to the overlap function 
jC{z), which is in general never trivial particularly when the number of terms is great. 
It is thus more popular to use non-squared from as in Eq. (|4-5p . 
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